import numpy as np
import pandas as pd
import math
from pyXSteam.XSteam import XSteam
import matplotlib.pyplot as plt
from matplotlib import rcParams
from tqdm import tqdm
from scipy.integrate import solve_ivp
from sklearn.metrics import mean_absolute_error, r2_score, root_mean_squared_error, mean_absolute_percentage_error

from sklearn.preprocessing import MinMaxScaler
from scipy.stats import gaussian_kde
import torch
import torch.nn as nn
from torch.utils.data import DataLoader, Dataset
from sklearn.model_selection import train_test_split
import time
import os

plt.rcParams['font.sans-serif'] = ['STZhongsong']
plt.rcParams['axes.unicode_minus'] = False

# =============================
# Configurations
# =============================
LOOKBACK = 30   # past 24 time steps as input
FORECAST_HORIZON = 10  # predict next 10 steps
CONFIDENCE = 0.9
TRAIN_SPLIT = 0.7
EPOCHS = 50
BATCH_SIZE = 32
LR = 0.001
DEVICE = 'cuda' if torch.cuda.is_available() else 'cpu'

SPEED_CHANGE_LIMIT = 16.5
SPEED_STEP = 1
TOLERANCE = 0.025
OUTPUT_FILE_OPT = "./results/冷端背压优化8月结果.xlsx"


# %%  汽轮机模型
# 汽轮机机理模型
class Steam_turbine:
    def __init__(self):
        self.Xsteam = XSteam(XSteam.UNIT_SYSTEM_MKS)  # 97 公式计算
        # 模型参数(初始值)
        self.K_list = [12.3 / 513.6, 6.9 / 583, 4 / 550, 2 / 520, 0.9 / 500, 0.53 / 485, 0.21 / 450, 0.08 / 420]  # 流量系数列表
        self.ita_c = 0.90  # 汽轮机末级效率  (0-1)
        self.ita_t = 0.90  # 汽轮机的机械效率 (0-1)

    def update_parameters(self, ita_c, ita_t, k_list):
        """更新参数"""
        self.ita_c = ita_c
        self.ita_t = ita_t
        self.K_list = k_list

    """汽轮机-调节级"""
    def regulating_level(self, Tst, Pst, Dst, k_list):
        """ 根据主蒸汽流量 Dst 的不同范围，计算调节级后的蒸汽参数和内能
        :param Tst: 主蒸汽温度
        :param Pst: 主蒸汽压力 (MPa)
        :param Dst: 主蒸汽流量 (t/h)
        :param k_list: 流量系数列表（用于计算压力）
        :return: Dt(调节级后蒸汽流量); pt(调节级后压力); ht(调节级后焓值); t(调节级后温度); ne(内能)
        """
        k_reg = k_list[0]  # 调节级流量系数
        # 常数
        pt_vwo = 14  # VWO（全阀开启）工况下调节级后压力 (MPa)
        pstb_vwo = 16.67  # VWO工况下部分开启调节阀后的压力 (MPa)
        Dstb_vwo = 290.76  # 设计工况下部分开启调节阀的蒸汽流量 (t/h)
        # pt_d = 12.3      # 设计工况下调节级后压力
        # Dst_d = 513.6    # 设计工况下主蒸汽流量
        fai = 1  # 喷嘴和动叶速度系数（无单位）
        Kn = 1.3  # 蒸汽绝热指数（无单位）
        Dst_vwo = 581.5138889  # VWO工况下主蒸汽流量 (t/h)
        hst = self.Xsteam.h_pt(Pst * 10, Tst)  # 主蒸汽焓值

        """调节级:根据主蒸汽流量 Dst 分为三种情况"""
        # """(低流量, 仅部分开启调节阀)"""
        if Dst <= 290.7569444:
            Dt = Dst  # 调节级后流量等于主蒸汽流量
            pt = Dst * k_reg  # 调节级后压力与流量成正比
            pst2 = Dst * ((pstb_vwo - pt_vwo) / Dstb_vwo) + pt  # 蒸汽流经部分开调节阀后的压力(部分开启阀门后的压力，通过线性插值计算)
            hst2 = self.Xsteam.h_pt(Pst * 10, Tst)  # 主蒸汽焓值
            vst2 = self.Xsteam.v_ph(pst2 * 10, hst2)  # 部分开启阀门后的比体积 (m³/kg)
            nm = Kn / (Kn - (Kn - 1) * (fai ** 2))  # 多变指数，考虑绝热膨胀
            vt2 = (pst2 / pt) ** (1 / nm) * vst2  # 通过部分开调节阀蒸汽流经喷嘴动叶后比体积
            tt2 = (1000 * pt * vt2 + 9.7 * pt - 123.1) / (0.461 + 0.01324 * pt)  # 通过部分开调节阀蒸汽流经喷嘴动叶后温度(经验公式计算温度)
            ht = self.Xsteam.h_pt(pt * 10, tt2)  # 调节级后焓值
            t = self.Xsteam.t_ph(pt * 10, ht)  # 调节级后温度

        # """(低流量, 仅部分开启调节阀)"""
        elif 290.7569444 < Dst <= 436.1354167:
            # 部分开启阀门
            pt = Dst * k_reg  # 变工况下调节级后压力
            pst2 = (Dst - 0.5 * Dst_vwo) * ((pstb_vwo - pt_vwo) / Dstb_vwo) + pt  # 蒸汽流经部分开调节阀后的压力
            hst2 = self.Xsteam.h_pt(Pst * 10, Tst)  # 主蒸汽焓值
            vst2 = self.Xsteam.v_ph(pst2 * 10, hst2)  # 部分开启阀门后的比体积 (m³/kg)
            nm = Kn / (Kn - (Kn - 1) * (fai ** 2))  # 多变指数，考虑绝热膨胀
            vt2 = (pst2 / pt) ** (1 / nm) * vst2  # 通过部分开调节阀蒸汽流经喷嘴动叶后比体积
            tt2 = (1000 * pt * vt2 + 9.7 * pt - 123.1) / (0.461 + 0.01324 * pt)  # 通过部分开调节阀蒸汽流经喷嘴动叶后温度
            ht2 = self.Xsteam.h_pt(pt * 10, tt2)  # 通过部分开启调节阀后蒸汽流经喷嘴动叶后的蒸汽焓值ht(部分开启阀门后焓值)
            # 全开阀门
            hst1 = self.Xsteam.h_pt(Pst * 10, Tst)  # 主蒸汽焓值
            pst1 = 0.95 * Pst  # 全开阀门后压力（假设 5% 压降）
            vst1 = self.Xsteam.v_ph(pst1 * 10, hst1)
            nm = Kn / (Kn - (Kn - 1) * (fai ** 2))  # 多变指数
            vt1 = (pst1 / pt) ** (1 / nm) * vst1  # 通过全开调节阀蒸汽流经喷嘴动叶后比体积
            tt1 = (1000 * pt * vt1 + 9.7 * pt - 123.1) / (0.461 + 0.01324 * pt)  # 通过全开调节阀蒸汽流经喷嘴动叶后温度
            ht1 = self.Xsteam.h_pt(pt * 10, tt1)  # 通过全开调节阀后蒸汽流经喷嘴动叶后的蒸汽焓值
            # 混合
            Dt = Dst  # 调节级后蒸汽流量
            ht = ((0.5 * Dst_vwo) * ht1 + (Dst - 0.5 * Dst_vwo) * ht2) / Dst  # 混合后调节级蒸汽焓值ht(按流量比例混合焓值)
            t = self.Xsteam.t_ph(pt * 10, ht)

            # 高流量(类似中流量，但混合比例调整为 75% 全开阀门流量)
        else:
            # 部分开启阀门
            pt = Dst * k_reg  # 变工况下调节级后压力
            pst2 = (Dst - 0.75 * Dst_vwo) * ((pstb_vwo - pt_vwo) / Dstb_vwo) + pt  # 蒸汽流经部分开调节阀后的压力
            hst2 = self.Xsteam.h_pt(Pst * 10, Tst)  # 主蒸汽焓值
            vst2 = self.Xsteam.v_ph(pst2 * 10, hst2)  # 流经部分开阀后的蒸汽的比体积
            nm = Kn / (Kn - (Kn - 1) * (fai ** 2))  # 多变指数，考虑绝热膨胀
            vt2 = (pst2 / pt) ** (1 / nm) * vst2  # 通过部分开调节阀蒸汽流经喷嘴动叶后比体积
            tt2 = (1000 * pt * vt2 + 9.7 * pt - 123.1) / (0.461 + 0.01324 * pt)  # 通过部分开调节阀蒸汽流经喷嘴动叶后温度
            ht2 = self.Xsteam.h_pt(pt * 10, tt2)  # 通过部分开启调节阀后蒸汽流经喷嘴动叶后的蒸汽焓值ht
            # 全开阀门
            hst1 = self.Xsteam.h_pt(Pst * 10, Tst)  # 主蒸汽焓值
            pst1 = 0.95 * Pst  # 蒸汽流经全开调节阀门后压力pst1
            vst1 = self.Xsteam.v_ph(pst1 * 10, hst1)
            nm = Kn / (Kn - (Kn - 1) * (fai ** 2))  # 多变指数
            vt1 = (pst1 / pt) ** (1 / nm) * vst1  # 通过全开调节阀蒸汽流经喷嘴动叶后比体积
            tt1 = (1000 * pt * vt1 + 9.7 * pt - 123.1) / (0.461 + 0.01324 * pt)  # 通过全开调节阀蒸汽流经喷嘴动叶后温度
            ht1 = self.Xsteam.h_pt(pt * 10, tt1)  # 通过全开调节阀后蒸汽流经喷嘴动叶后的蒸汽焓值
            # 混合
            Dt = Dst  # 调节级后蒸汽流量
            ht = ((0.75 * Dst_vwo) * ht1 + (Dst - 0.75 * Dst_vwo) * ht2) / Dst  # 混合后调节级蒸汽焓值ht
            t = self.Xsteam.t_ph(pt * 10, ht)
        # 内能 (kW)，表示调节级能量转换。
        ne = Dst * (hst - ht)
        return Dt, pt, ht, t, ne

    """汽轮机-压力级"""
    "压力级第一级"

    def pressure_level_1(self, Dt, pt, ht, k_list):
        """
        :param Dt: 入口蒸汽流量 (t/h)
        :param pt: 入口压力 (MPa)
        :param ht: 入口焓值 (kJ/kg)
        :param k_list: 流量系数列表
        :return: D1(抽汽流量); p1(抽汽压力); h1(抽汽焓值); ne(内能)
        """
        k_pre_1 = k_list[1]
        p1d = 6.90877688  # 设计工况抽汽压力
        xl = 0.9  # 级间效率
        Dj1d = 583.3333333  # 设计工况下第一级蒸汽流量
        p1 = Dt * k_pre_1  # 第一级抽汽压力，即第一级出口蒸汽压力(抽汽压力与流量成正比)
        s1 = self.Xsteam.s_ph(pt * 10, ht)  # 级组入口蒸汽熵
        h11 = self.Xsteam.h_ps(p1 * 10, s1)  # 第一级抽汽理想焓值(等熵膨胀)
        h1 = ht - xl * (ht - h11)  # 第一级实际抽汽焓值(实际焓值，考虑效率)
        Pa = 0.101  # 大气压 (MPa)
        EGJCQ = 3
        ps1 = (p1 + Pa) * (1 - EGJCQ / 100) - Pa  # 抽汽口压力，考虑 3% 损失
        if ps1 > p1:
            D1 = 0
        else:
            D1 = math.sqrt(p1 - ps1) * (-0.4972 * p1 ** 3 + 7.227 * p1 ** 2 - 19.44 * p1 + 48.69)  # 第一级抽汽流量
        ne = Dt * (ht - h1)  # 内能
        return D1, p1, h1, ne

    "压力级第二级"

    def pressure_level_2(self, D1, p1, h1, Dt, k_list):
        k_pre_2 = k_list[2]
        p2d = 4.065  # 设计工况抽汽压力
        xl = 0.9  # 级间效率
        Dj2d = 550.548348  # 设计工况下第二级蒸汽流量
        Dj2 = Dt - D1  # 第二级蒸汽流量
        p2 = Dj2 * k_pre_2  # 第二级抽汽压力，即第二级出口蒸汽压力
        s2 = self.Xsteam.s_ph(p1 * 10, h1)  # 级组二入口蒸汽熵
        h22 = self.Xsteam.h_ps(p2 * 10, s2)  # 第二级抽汽理想焓值
        h2 = h1 - xl * (h1 - h22)  # 第二级实际抽汽焓值
        Pa = 0.101
        EGJCQ = 3
        ps2 = (p2 + Pa) * (1 - EGJCQ / 100) - Pa  # 抽汽口压力，考虑 3% 损失
        if ps2 > p2:
            D2 = 0
        else:
            D2 = math.sqrt(p2 - ps2) * (-2.75 * p2 ** 3 + 23.6 * p2 ** 2 - 38.75 * p2 + 62.27)
        ne = Dj2 * (h1 - h2)
        return D2, p2, h2, Dj2, ne

    "压力级第三级"

    def pressure_level_3(self, Pre, Tre, Dre, k_list):
        """
        :param Pre: 再热蒸汽压力
        :param Tre: 再热蒸汽温度
        :param Dre: 再热蒸汽流量
        :param k_list: 流量系数
        :return:
        """
        k_pre_3 = k_list[3]
        p3d = 2.090  # 设计工况下第三级抽汽压力
        xl = 0.9  # 级间效率
        Dj3d = 520.07  # 设计工况下第三级蒸汽流量
        p3 = Dre * k_pre_3  # 第三级抽汽压力，即第三级出口蒸汽压力
        hre = self.Xsteam.h_pt(Pre * 10, Tre)  # 级组三入口蒸汽焓值
        s3 = self.Xsteam.s_ph(Pre * 10, hre)  # 级组三入口蒸汽熵值
        h33 = self.Xsteam.h_ps(p3 * 10, s3)  # 第三级抽汽理想焓值
        h3 = hre - xl * (hre - h33)  # 第三级实际抽汽焓值
        Pa = 0.101
        EGJCQ = 3
        ps3 = (p3 + Pa) * (1 - EGJCQ / 100) - Pa  # 抽汽口压力，考虑 3% 损失
        if ps3 > p3:
            D3 = 0
        else:
            D3 = math.sqrt(p3 - ps3) * (-6.841 * p3 ** 3 + 31.39 * p3 ** 2 - 17.2 * p3 + 22.72)  # 第三级抽汽流量
        Dj3 = Dre - D3  # 第三级蒸汽流量
        ne = Dre * (hre - h3)
        return D3, p3, h3, Dj3, ne, hre

    "压力级第四级"

    def pressure_level_4(self, h3, p3, Dj3, k_list):
        k_pre_4 = k_list[4]
        p4d = 0.916990  # 设计工况下第四级抽汽压力
        xl = 0.9  # 级间效率
        Dj4d = 500.07  # 设计工况下第四级蒸汽流量
        p4 = Dj3 * k_pre_4  # 第四级抽汽压力，即第四级出口蒸汽压力
        s4 = self.Xsteam.s_ph(p3 * 10, h3)  # 级组四入口蒸汽熵
        h44 = self.Xsteam.h_ps(p4 * 10, s4)  # 第四级抽汽理想焓值
        h4 = h3 - xl * (h3 - h44)  # 第四级实际抽汽焓值
        Pa = 0.101
        EGJCQ = 3
        ps4 = (p4 + Pa) * (1 - EGJCQ / 100) - Pa  # 抽汽口压力，考虑 3% 损失
        if ps4 > p4:
            D4 = 0
        else:
            D4 = math.sqrt(p4 - ps4) * (-64.69 * p4 ** 3 + 144.7 * p4 ** 2 - 47.91 * p4 + 35.01)  # 第四级抽汽流量
        ne = Dj3 * (h3 - h4)
        Dj4 = Dj3 - D4  # 第四级蒸汽流量

        return D4, p4, h4, Dj4, ne

    "压力级第五级"

    def pressure_level_5(self, h4, p4, Dj4, k_list):
        k_pre_5 = k_list[5]
        p5d = 0.53591  # 设计工况下第五级抽汽压力
        xl = 0.9  # 级间效率
        Dj5d = 485.54  # 设计工况下第五级蒸汽流量
        p5 = Dj4 * k_pre_5  # 第五级抽汽压力，即第五级出口蒸汽压力
        s5 = self.Xsteam.s_ph(p4 * 10, h4)  # 级组五入口蒸汽熵
        h55 = self.Xsteam.h_ps(p5 * 10, s5)  # 第五级抽汽理想焓值
        h5 = h4 - xl * (h4 - h55)  # 第五级实际抽汽焓值
        Pa = 0.101
        EGJCQ = 3
        ps5 = (p5 + Pa) * (1 - EGJCQ / 100) - Pa  # 抽汽口压力，考虑 3% 损失
        if ps5 > p5:
            D5 = 0
        else:
            D5 = math.sqrt(p5 - ps5) * (-365.1 * p5 ** 3 + 529.1 * p5 ** 2 - 68.05 * p5 + 56.83)  # 第五级抽汽流量
        ne = Dj4 * (h4 - h5)
        Dj5 = Dj4 - D5  # 第五级蒸汽流量

        return D5, p5, h5, Dj5, ne

    "压力级第六级"

    def pressure_level_6(self, h5, p5, Dj5, k_list):
        k_pre_6 = k_list[6]
        p6d = 0.21902712  # 设计工况下第六级抽汽压力
        xl = 0.9  # 级间效率
        Dj6d = 450.07  # 设计工况下第六级蒸汽流量
        p6 = Dj5 * k_pre_6  # 第六级抽汽压力，即第六级出口蒸汽压力
        s6 = self.Xsteam.s_ph(p5 * 10, h5)  # 级组六入口蒸汽熵
        h66 = self.Xsteam.h_ps(p6 * 10, s6)  # 第六级抽汽理想焓值
        h6 = h5 - xl * (h5 - h66) + 110  # 第六级实际抽汽焓值
        Pa = 0.101
        EGJCQ = 3
        ps6 = (p6 + Pa) * (1 - EGJCQ / 100) - Pa  # 抽汽口压力，考虑 3% 损失
        if ps6 > p6:
            D6 = 0
        else:
            D6 = math.sqrt(p6 - ps6) * (9233 * p6 ** 3 - 4349 * p6 ** 2 + 1313 * p6 + 2.72)  # 第六级抽汽流量
        ne = Dj5 * (h5 - h6)
        Dj6 = Dj5 - D6  # 第六级蒸汽流量
        return D6, p6, h6, Dj6, ne

    "压力级第七级"

    def pressure_level_7(self, h6, p6, Dj6, k_list):
        k_pre_7 = k_list[7]
        p7d = 0.081204  # 设计工况下第七级抽汽压力
        xl = 0.9  # 级间效率
        Dj7d = 420.95  # 设计工况下第七级蒸汽流量
        p7 = Dj6 * k_pre_7  # 第七级抽汽压力，即第七级出口蒸汽压力
        s7 = self.Xsteam.s_ph(p6 * 10, h6)  # 级组七入口蒸汽熵
        h77 = self.Xsteam.h_ps(p7 * 10, s7)  # 第七级抽汽理想焓值
        h7 = h6 - xl * (h6 - h77)  # 第七级实际抽汽焓值
        Pa = 0.101
        EGJCQ = 3
        ps7 = (p7 + Pa) * (1 - EGJCQ / 100) - Pa  # 抽汽口压力，考虑 3% 损失
        if ps7 > p7:
            D7 = 0
        else:
            D7 = math.sqrt(p7 - ps7) * (5608 * p7 ** 3 - 604 * p7 ** 2 + 2841 * p7 - 54.62)  # 第七级抽汽流量
        ne = Dj6 * (h6 - h7)
        Dj7 = Dj6 - D7  # 第七级蒸汽流量
        return D7, p7, h7, Dj7, ne

    "末级"

    def final_level(self, h7, p7, Dj7, Pc):
        """
        :param h7: 压力级出口焓值（第七级焓值）
        :param p7: 压力级出口压力（第七级压力）
        :param Dj7: 压力级流量（第七级流量）
        :param Pc: 背压
        :return: hc（排汽焓值）; Dc (排汽流量); ne(内能); Tp(排汽温度); Tl_p(饱和温度)
        """
        # 汽轮机末级效率
        xl = self.ita_c
        Dc = Dj7  # 末级蒸汽流量
        sm = self.Xsteam.s_ph(p7 * 10, h7)  # 末级入口蒸汽熵
        hm = self.Xsteam.h_ps(Pc * 10, sm)  # 末级抽汽理想焓值
        hc = h7 - xl * (h7 - hm)  # 末级排汽焓值
        ne = Dc * (h7 - hc)  # 内能
        Tc = self.Xsteam.t_ph(Pc * 10, hc)  # 排汽温度

        return hc, Dc, Tc, ne

    "汽轮机仿真：综合模拟汽轮机各阶段，计算负荷"

    def Turbine_simulation(self, Tst, Pst, Dst, Pre, Tre, Pc):
        """
        :param Tst: 主蒸汽温度
        :param Pst: 主蒸汽压力
        :param Dst: 主蒸汽流量
        :param Pre: 再热蒸汽压力
        :param Tre: 再热蒸汽温度
        :param Pc: 背压
        :return: Dt(调节级后流量); Pt(调节级压力); Dc(末级排汽流量); Ne(总负荷)
        """
        k_list = self.K_list  # 流量系数列表
        # 调节级
        Dt, Pt, ht, t, ne_regulate = self.regulating_level(Tst, Pst, Dst, k_list)
        # 压力级第一级
        D1, p1, h1, ne1 = self.pressure_level_1(Dt, Pt, ht, k_list)
        # 压力级第二级
        D2, p2, h2, Dj2, ne2 = self.pressure_level_2(D1, p1, h1, Dt, k_list)
        # 压力级第三级 (再热器)
        D3, p3, h3, Dj3, ne3, hre = self.pressure_level_3(Pre, Tre, Dj2, k_list)
        # 压力级第四级
        D4, p4, h4, Dj4, ne4 = self.pressure_level_4(h3, p3, Dj3, k_list)
        # 压力级第五级
        D5, p5, h5, Dj5, ne5 = self.pressure_level_5(h4, p4, Dj4, k_list)
        # 压力级第六级
        D6, p6, h6, Dj6, ne6 = self.pressure_level_6(h5, p5, Dj5, k_list)
        # 压力级第七级
        D7, p7, h7, Dj7, ne7 = self.pressure_level_7(h6, p6, Dj6, k_list)
        # 末级
        hc, Dc, Tc, ne_final = self.final_level(h7, p7, Dj7, Pc)
        # 总负荷 (总内能乘以机械效率并转换为 MW)
        Ne = (ne_regulate + ne1 + ne2 + ne3 + ne4 + ne5 + ne6 + ne7 + ne_final) * self.ita_t
        Ne = Ne * 0.001  # 转换为 MW

        return Dt, Pt, Dc, hc, Tc, Ne


# %% 直接空冷凝汽器的机理模型
class Air_Condenser:
    def __init__(self):
        self.Xsteam = XSteam(XSteam.UNIT_SYSTEM_MKS)
        # 凝汽器的辨识参数
        self.kc = 0.04  # 凝汽器换热系数
        self.cf = 0.9  # 风机集群因子

        # 凝汽器的模型参数
        self.Vn = 2.5  # 额定风速 (m/s)
        self.Ac = 1440216  # 换热表面积 (m²)
        self.Af = 10731  # 风机覆盖面积 (m²)
        self.rho_air = 1.043  # 空气密度 (kg/m³)
        self.Cp_air = 1.005  # 空气比热容 (kJ/(kg·℃))
        self.Qvd = 550  # 风机设计体积流量 (m³/s)
        self.n_design = 110  # 设计转速
        self.Rs = 461.5  # 水蒸气气体常数 (J/(kg·K))
        self.Vc = 9600 * 56  # 凝汽器体积 (m³)
        self.kec = 1.5  # 凝结水流量修正系数

    def update_parameters(self, kc, cf):
        """更新参数"""
        self.kc = kc
        self.cf = cf

    def Air_Condenser_Sim(self, t, x, u):
        """定义空冷凝汽器的微分方程，用于 solve_ivp 求解"""
        Pc = x[0]  # 背压 (MPa)
        Dc, hc, Ta, fan_speed, fan_number = u  # 控制输入：排汽流量、排汽焓值、环境温度、风机转速、风机运行台数

        k = self.kc
        cluster_factor = self.cf

        # 凝汽器的模型参数
        Vn = self.Vn
        Ac = self.Ac
        Af = self.Af
        rho_air = self.rho_air
        Cp_air = self.Cp_air
        Qvd = self.Qvd
        n_design = self.n_design
        Rs = self.Rs
        Vc = self.Vc
        kec = self.kec

        pc = Pc * 10  # MPa to bar
        hc_sa = self.Xsteam.hV_p(pc)  # 凝汽器压力对应的饱和蒸汽焓(理想排汽焓值)
        Tc_sa = self.Xsteam.tsat_p(pc)  # 凝汽器压力对应的饱和温度
        hlc = self.Xsteam.hL_p(pc)  # 凝汽器压力对应的饱和水焓

        # 计算风机转速
        Qv = Qvd * (fan_speed / n_design)  # 单台风机体积流量
        Vf_current = (Qv * fan_number * cluster_factor) / Af  # 当前风速
        Vf = Vf_current

        kc = k * (np.power(complex(Vf / Vn), 0.6)).real  # 修正后的传热系数
        NTU = (kc * Ac) / (rho_air * Vf * Af * Cp_air)  # 传热单元数NTU
        Qc = (rho_air * Vf * Af * Cp_air) * (1 - math.exp(-NTU)) * (Tc_sa - Ta)  # 此汽轮机排汽和环境空气之间的换热量 Qc
        Dec = Qc / (hc - hlc) * kec  # 凝汽器内凝结的排汽质量流量
        dPc_dt = (Rs * (Tc_sa + 273.15) * (Dc - Dec)) / Vc  # Pc 的更新量（单位为：Pa）
        dPc_dt = dPc_dt * 0.000001  # Pa → MPa

        return np.array([dPc_dt])

    def Calculate_current_state(self, Pc, Dc, hc, Ta, fan_speed, fan_number):
        """计算当前时刻的凝结水流量和风速"""
        pc = Pc * 10  # MPa to bar
        hc_sa = self.Xsteam.hV_p(pc)  # 凝汽器压力对应的饱和蒸汽焓(理想排汽焓值)
        Tc_sa = self.Xsteam.tsat_p(pc)  # 凝汽器压力对应的饱和温度
        hlc = self.Xsteam.hL_p(pc)  # 凝汽器压力对应的饱和水焓

        k = self.kc
        cluster_factor = self.cf
        Vn = self.Vn
        Ac = self.Ac
        Af = self.Af
        rho_air = self.rho_air
        Cp_air = self.Cp_air
        Qvd = self.Qvd
        n_design = self.n_design
        kec = self.kec

        # 计算风机转速
        Qv = Qvd * (fan_speed / n_design)  # 单台风机体积流量
        Vf_current = (Qv * fan_number * cluster_factor) / Af  # 当前风速
        Vf = Vf_current

        kc = k * (np.power(complex(Vf / Vn), 0.6)).real  # 修正后的传热系数
        NTU = (kc * Ac) / (rho_air * Vf * Af * Cp_air)  # 传热单元数NTU
        Qc = (rho_air * Vf * Af * Cp_air) * (1 - math.exp(-NTU)) * (Tc_sa - Ta)  # 此汽轮机排汽和环境空气之间的换热量 Qc
        Dec = Qc / (hc - hlc) * kec  # 凝汽器内凝结的排汽质量流量

        return Vf, Dec, Pc, kc

    def Condenser_Simulation(self, Pc, Dc, hc, Ta, fan_speed, fan_number, t_end=60.0, dt=0.1):
        """ 执行一步积分仿真，返回当前时刻参数以及积分后下一时刻的背压
            t_end: 仿真时间段长度 (默认60秒，即采样时间)
            dt: 积分步长 (默认0.1秒)
        """
        # 计算当前时刻的状态(风速；凝结水流量；背压；kc参数)
        Vf, Dec, Pc, kc = self.Calculate_current_state(Pc, Dc, hc, Ta, fan_speed, fan_number)

        tspan = np.arange(0, t_end + dt, dt)
        u_ctrl = np.array([Dc, hc, Ta, fan_speed, fan_number])
        state = np.array([Pc])

        # 背压模型进行积分
        sol = solve_ivp(
            fun=lambda t, x: self.Air_Condenser_Sim(t, x, u_ctrl),
            t_span=(tspan[0], tspan[-1]),
            y0=state,
            t_eval=tspan)
        Pc_next = sol.y[:, -1]  # 取最后一个作为下一时刻的状态
        Pc_next = Pc_next[0]

        return Pc_next, Vf, Dec, kc, Pc


# %% 风机功耗模型
class Air_Power:
    def __init__(self):
        pass

    @staticmethod
    def Air_Power_Sim(fan_speed, fan_number):
        # 计算风机功耗
        n0 = 110  # 额定风机转速
        P0 = 115  # 额定风机功率(115KW)
        Power_fan = ((fan_speed / n0) ** 3) * P0 * fan_number * 0.001  # 单位:MW
        return Power_fan


# %% 耦合仿真
class Model_Simulation:
    def __init__(self, Datas, speed_change_limit, speed_step, output_filepath, tolerance=0.025):
        """ 使用数据和参数初始化模拟类，汽轮机仿真 """
        self.data = Datas  # 全部待仿真数据
        self.N = len(self.data)-11  # 数据长度，留十个作为优化的未来
        self.tolerance = tolerance  # 稳态判断容忍度

        # 真实数据
        self.Ne_data = self.data[:, 0]              # 机组负荷
        self.Pc_data = self.data[:, 1]              # 背压 (单位：kPa)
        self.Pst_data = self.data[:, 2]             # 主汽压力
        self.Dst_data = self.data[:, 3]             # 主汽流量
        self.Tst_data = self.data[:, 4]             # 主汽温度
        self.Feed_water_data = self.data[:, 5]      # 主给水流量
        self.Pre_data = self.data[:, 6]             # 再热蒸汽压力
        self.Tre_data = self.data[:, 7]             # 再热蒸汽温度
        self.Pt_data = self.data[:, 14]             # 调节级压力
        self.Tc_data = self.data[:, 13]             # 末级排汽温度
        self.P_pre_data = self.data[:, -7:]         # 各个压力级的抽汽压力
        self.Ta_data = self.data[:, 10]             # 环境温度
        self.fan_speed_data = self.data[:, 11]      # 风机转速
        self.fan_number_data = self.data[:, 12]     # 风机台数
        self.Dec_data = self.data[:, 9]             # 凝结水流量
        self.fan_power_data = self.data[:, 15]      # 风机功耗
        self.net_power_data = self.data[:, 16]      # 净发电功率

        self.Xsteam = XSteam(XSteam.UNIT_SYSTEM_MKS)
        self.Turbine = Steam_turbine()                  # 蒸汽轮机模型
        self.Condenser = Air_Condenser()                # 空冷凝汽器模型
        self.Condenser_Power = Air_Power()              # 风机功耗模型

        # 结果保存
        self.Pt_sim = []            # 调节级压力仿真
        self.Ne_sim = []            # 机组负荷仿真
        self.Dc_sim = []            # 末级排汽流量仿真
        self.Tc_sim = []            # 末级排汽温度仿真
        self.hc_sim = []            # 末级排汽焓值仿真
        self.Pc_sim = []            # 背压仿真
        self.Dec_sim = []           # 凝结水流量仿真

        self.ita_t_sim = []
        self.ita_c_sim = []
        self.kc_sim = []
        self.cf_sim = []
        self.Rf_sim = []

        # 在线梯度下降更新效率参数，滑动窗口最小二乘法更新流量系数
        self.ita_c = self.Turbine.ita_c
        self.ita_t = self.Turbine.ita_t
        self.k_list = self.Turbine.K_list
        self.learning_rate = 2.5e-6  # 学习率
        self.P_buf = []
        self.D_buf = []

        self.kc = self.Condenser.kc
        self.cf = self.Condenser.cf

        self.Pc_current = None  # 初始化当前时刻背压

        # 优化相关
        self.prev_optimal_speed=None
        self.speed_change_limit=speed_change_limit
        self.speed_step=speed_step
        self.current_sample_results = []  # 存储当前样本在所有风机转速下的结果
        self.all_samples_results = []  # 所有样本在所有转速下的结果
        self.result_list = [] # 存储所有样本的最优结果
        self.output_filepath = output_filepath

        self.load_error=0
        self.pc_error=0
        self.pt_error=0
        self.dec_error=0

        self.opt_net_power=0
        self.opt_ne = 0
        self.opt_pc = 0
        self.opt_fan_speed = 0
        self.opt_fan_power = 0
        self.opt_dst = 0
        self.opt_coal = 0
        self.opt_coal_cost = 0
        self.opt_fan_cost =0
        self.opt_profit=0

    def Sim_Step(self, input_data, Pc):
        """ 执行单步模拟仿真 """
        Tst, Pst, Dst, Pre, Tre, Ta, fan_speed, fan_number = input_data

        Dt, Pt, Dc, hc, Tc, Ne = self.Turbine.Turbine_simulation(Tst, Pst, Dst, Pre, Tre, Pc)
        Pc_next, Vf, Dec, kc, Pc = self.Condenser.Condenser_Simulation(Pc, Dc, hc, Ta, fan_speed, fan_number)

        return np.array([Dt, Pt, Dc, hc, Tc, Ne, Pc, Vf, Dec, Pc_next])

    @staticmethod
    def Net_power(Ne, n_opt, n_num):
        power_consumption_fan = (n_opt / 110) ** 3 * 115 * 0.001 * n_num
        current_net_power = Ne - power_consumption_fan
        return power_consumption_fan, current_net_power

    def optimize_for_speed(self, mech_in, m1, agc, Pc_k, fan_power, k_Dt, k_coal, tolerance, k=None):
        """
            优化风机转速，根据仿真负荷与实际负荷的差值动态调整主汽流量

            参数:
            agc: self.Ne_data[k]
            Pc_k: self.Pc_data[k]
            fan_power: self.fan_power_data[k]
            k_Dt = 0.6039  # 主汽流量与给煤量关系参数，需从外部输入
            k_coal = 0.9574  # 负荷与主汽流量关系参数，需从外部输入
            max_iterations: 负荷-主汽流量迭代的最大次数
            """
        # k_coal = 0.9574  # 负荷与主汽流量关系参数，需从外部输入
        # b = -20.349  # 负荷与主汽流量关系参数，需从外部输入
        # k_Dt = 0.6039  # 主汽流量与给煤量关系参数，需从外部输入
        # d = 29.2676  # 主汽流量与给煤量关系参数，需从外部输入

        # 输入煤价和电价
        # load_rate = data1[1]/600  # 实际负荷
        # coal_price = 500 if load_rate <= 0.5 else 550 if load_rate <= 0.75 else 625
        coal_price = 550.0 / 1000  # 煤价 (元/吨)
        electricity_price = 0.35291  # 电价 (元/kWh)
        max_iterations = 50

        # 保存原始数据
        # actual_load = data1[1]  # 实际负荷 self.Ne_data[k]
        # actual_back_pressure = data1[2]  # kPa，实际背压 self.Pc_data[k] self.Pc_current
        actual_D = mech_in[2]  # 原始主汽流量
        actual_fan_speed = mech_in[-2]  # 实际风机转速
        # 当前输入
        mech_in[-2] = m1
        n_num = mech_in[-1]  # 风机运行台数

        # 主汽流量-负荷迭代调整
        for iteration in range(max_iterations):
            # 运行仿真计算当前背压和负荷
            # 单步仿真
            states = self.Sim_Step(mech_in, Pc_k)
            Pc_k = states[-1]
            current_Pc = states[-1] * 1000  # 当前背压 (MPa-→kPa)
            current_Ne = states[5]  # 当前负荷 (MW)

            delta_Ne = agc - current_Ne  # 仿真与实际负荷差
            print(f"agc: {agc}, current_Ne: {current_Ne}, delta_Ne: {delta_Ne}")
            delta_Ne = np.clip(delta_Ne, -20, 20)

            dD = (k_Dt * delta_Ne)

            # 根据负荷差值调整主汽流量
            D_adj = mech_in[2] + dD

            # 更新数据中的主汽流量
            mech_in[2] = D_adj

            # 判断是否收敛
            if abs(delta_Ne / agc) < tolerance or iteration == max_iterations - 1:
                # print(f"风机转速 {m1}: 负荷={current_Ne:.2f} MW, 背压={current_Pc:.2f} kPa, "
                #       f"主汽流量={D_adj:.2f} t/h, 负荷差值={load_diff:.2f} MW")# 打印迭代信息
                break

        # 最终背压和负荷
        final_back_pressure = current_Pc
        final_load = current_Ne

        # 基于主汽流量变化计算给煤量减少量，给煤量与主汽流量成线性关系: Coal = c * D + d
        coal_reduction = k_coal * (actual_D - D_adj)  # 给煤量减少量
        print(f"actual_D: {actual_D}, D_adj: {D_adj}")
        print(f"coal reduction: {coal_reduction} k_coal: {k_coal}")

        # 计算给煤节省的成本 (元/分钟 或 元)
        coal_cost_saving = coal_reduction * 60 * coal_price
        print(f"coal cost saving: {coal_cost_saving}")

        # 计算净功率
        power_consumption_fan, net_power = self.Net_power(float(final_load), m1, n_num)
        print(f"power_consumption_fan: {power_consumption_fan}")

        # 计算风机增加的功耗成本 (元/分钟)
        power_cost_increase = (power_consumption_fan - fan_power) * 1000 / 60 * electricity_price  # 转换为 kW

        print(f"fan_power: {fan_power}")
        print(f"{m1}:coal_cost_saving-{coal_cost_saving}, power_cost_increase-{power_cost_increase}")
        # 经济效益计算
        objective_value = coal_cost_saving - power_cost_increase

        if k is None:
            return ([
                round(net_power, 2),  # 最大净功率
                round(final_load, 2),  # 优化负荷
                round(final_back_pressure, 2),  # 最优背压
                round(m1, 2),  # 优化风机转速
                round(power_consumption_fan, 2),  # 优化风机功耗
                round(D_adj, 2),  # 调整后的主汽流量
                round(coal_reduction, 2),  # 给煤量减少量
                round(coal_cost_saving, 2),  # 给煤节省成本 (元/分钟)
                round(power_cost_increase, 2),  # 风机增加成本 (元/分钟)
                round(objective_value, 2)  # 目标函数值 (元/分钟)
            ],states)
        else:
            return ([
                    round(self.Ta_data[k], 2),  # 环境温度
                    round(net_power, 2),  # 最大净功率
                    round(self.net_power_data[k+1], 2),  # 实际净功率
                    round(final_load, 2),  # 优化负荷
                    round(self.Ne_data[k+1], 2),  # 实际负荷
                    round(final_back_pressure, 2),  # 优化背压
                    round(self.Pc_data[k+1], 2),  # 实际背压
                    round(m1, 2),  # 优化风机转速
                    round(self.fan_speed_data[k+1], 2),  # 实际风机转速
                    round(self.fan_number_data[k], 2),  # 风机运行台数
                    round(power_consumption_fan, 2),  # 优化风机功耗
                    round(self.fan_power_data[k+1], 2),  # 实际风机功耗
                    round(D_adj, 2),  # 调整后的主汽流量
                    round(self.Dst_data[k+1], 2),  # 原始主汽流量
                    round(coal_reduction, 2),  # 给煤量减少量
                    round(coal_cost_saving, 2),  # 给煤节省成本 (元/分钟)
                    round(power_cost_increase, 2),  # 风机增加成本 (元/分钟)
                    round(objective_value, 2)  # 目标函数值 (元/分钟)
                ]
            , states)

    def Run_sim_opt(self):
        """ 仿真运行 然后 优化 """
        self.Pc_current = self.Pc_data[0] / 1000  # 初始状态 (Mpa)
        Pc_pre = self.Pc_current

        for k in tqdm(range(self.N)):
            # 当前输入
            mech_in = [self.Tst_data[k], self.Pst_data[k], self.Dst_data[k], self.Pre_data[k], self.Tre_data[k],
                       self.Ta_data[k], self.fan_speed_data[k], self.fan_number_data[k]]

            # 使用当前参数运行模型
            self.Turbine.update_parameters(self.ita_c, self.ita_t, self.k_list)
            self.Condenser.update_parameters(self.kc, self.cf)

            # 单步仿真
            states = self.Sim_Step(mech_in, self.Pc_current)

            # 保存仿真结果
            self.Pt_sim.append(states[1])
            self.Dc_sim.append(states[2])
            self.hc_sim.append(states[3])
            self.Tc_sim.append(states[4])
            self.Ne_sim.append(states[5])
            self.Pc_sim.append(states[6] * 1000)
            self.Dec_sim.append(states[8])
            Pc_next = states[-1]

            # 更新模型参数(在线梯度下降方法)
            p_real = np.hstack([self.Pt_data[k], self.P_pre_data[k, :]])
            self.Update_Turbine_Parameters(mech_in, p_real, self.Ne_data[k], self.Pc_current)

            if k > 0:
                mech_pre = [self.Tst_data[k - 1], self.Pst_data[k - 1], self.Dst_data[k - 1], self.Pre_data[k - 1], self.Tre_data[k - 1],
                            self.Ta_data[k - 1], self.fan_speed_data[k - 1], self.fan_number_data[k - 1]]
                self.Update_Condenser_Parameters(mech_pre, Pc_pre, self.Pc_data[k])

            self.ita_c_sim.append(self.ita_c)
            self.ita_t_sim.append(self.ita_t)
            self.kc_sim.append(self.kc * 1000)
            self.Rf_sim.append(1 / (self.kc * 1000) - 1 / 40.5) #清洁因子？
            self.cf_sim.append(self.cf)

            Pc_pre = self.Pc_current
            self.Pc_current = Pc_next  # Pc用仿真值一直迭代
            # TODO: optimize alg
            start_t = time.time()

            # ---- 滑动窗口线性回归辨识k_coal和k_Dt
            if k > 360:
                # 17是给煤量
                k_coal = (self.data[k - 360:k, 17].reshape(1, -1) @ np.linalg.pinv(
                    np.concatenate((self.data[k - 360:k, 3].reshape(1, -1), np.ones((1, 360))))))  # .reshape(-1)
                k_Dt = (self.data[k - 360:k, 3].reshape(1, -1) @ np.linalg.pinv(
                    np.concatenate((self.data[k - 360:k, 0].reshape(1, -1), np.ones((1, 360))))))  # .reshape(-1)
                k_coal = k_coal[0, 0]
                k_Dt = k_Dt[0, 0]
            else:
                k_coal = 0.1929
                k_Dt = 0.8624

            # ---- 根据负荷率计算背压上下限 ----
            full_load = 600
            load_rate = self.Ne_data[k] / full_load  # 当前数据点的负荷率

            max_back_pressure = 41.8182 * load_rate + 16.6364 if 0.2 < load_rate <= 0.75 else 25 if load_rate <= 0.2 else 48
            min_back_pressure = 8.75 * load_rate - 0.75 if load_rate > 0.1 else 0
            # min_speed, max_speed = 44, 121

            # 确定风机转速搜索范围
            if self.prev_optimal_speed is None:
                # 第一个样本使用完整范围
                min_speed, max_speed = 44, 120.5
                # min_speed, max_speed = 70, 100
            else:
                # 后续样本考虑相邻时刻变化率约束
                min_speed = max(44, self.prev_optimal_speed - self.speed_change_limit)
                max_speed = min(120.5, self.prev_optimal_speed + self.speed_change_limit)

            print(f"样本 {k + 1}/{self.N}: 风机转速搜索范围 [{min_speed:.1f}, {max_speed:.1f}] rpm")

            # 生成要遍历的风机转速列表
            speed_range = np.arange(min_speed, max_speed + self.speed_step, self.speed_step)
            print(f"需要评估的风机转速点数量: {len(speed_range)}")

            # 存储当前样本在所有风机转速下的结果
            #current_sample_results = []
            best_result = None
            best_objective = -float('inf')
            best_speed = None

            # 遍历每个风机转速，评估经济效益
            prev_valid_bp = float('inf')  # 初始化上一个有效背压为无穷大
            for m1 in speed_range:
                data_copy = mech_in.copy()
                result, state_result = self.optimize_for_speed(mech_in=data_copy, m1=m1, agc=self.Ne_data[k], Pc_k=self.Pc_data[k]/1000, fan_power=self.fan_power_data[k], k_Dt=k_Dt, k_coal=k_coal, tolerance=self.tolerance, k=k)
                if m1 > min_speed:
                    bp = result[5]  # 最优背压 (kPa)
                    objective = result[-1]  # 目标函数值 (元/分钟)

                    # 添加风机转速和样本索引到结果
                    extended_result = [k + 1] + result  # 添加样本编号和风机转速
                    self.current_sample_results.append(extended_result)

                    # 检查背压是否在合理范围内
                    if min_back_pressure <= bp <= max_back_pressure:
                        # 检查背压是否不再下降
                        # if bp >= prev_valid_bp:
                        #     print(f"风机转速为{m1:.2f} 时，背压{bp:.2f} kPa不再下降（上一个有效背压: {prev_valid_bp:.2f} kPa）")
                        #     break  # 背压不再下降时跳出循环

                        # prev_valid_bp = bp  # 更新上一个有效背压

                        if objective > best_objective:
                            best_objective = objective
                            best_result = extended_result
                            best_speed = m1
                    else:
                        # print(f"风机转速为{m1:.2f} 时，背压{bp:.2f} kPa超出合理范围")
                        continue  # 背压超出范围时跳出循环
                    print(f"{m1}:bp-{bp}, obj-{objective}")

                # 打印进度
                if (m1 - min_speed) % (self.speed_step * 10) == 0:
                    print(f"  已评估 {int((m1 - min_speed) / self.speed_step + 1)}/{len(speed_range)} 个转速点")

            # 将当前样本的所有结果添加到总结果列表
            self.all_samples_results.extend(self.current_sample_results)

            end_t = time.time()
            print(f'样本 {k + 1}/{self.N} 优化完成，耗时{(end_t - start_t):.2f}秒，'
                  f'最优风机转速: {best_speed} rpm，经济效益: {best_objective:.2f} 元/分钟')

            # 记录最优结果
            if best_result is not None:
                print(f"已找到最优解（数据点 {k + 1}/{self.N}）")

                self.result_list.append(best_result)
                self.prev_optimal_speed = best_speed  # 保存最优风机转速

            if k % 360 == 0 or k == self.N - 1:

                # 保存所有样本在所有风机转速下的完整结果
                all_results_df = pd.DataFrame(self.all_samples_results, columns=[
                    '样本编号',
                    '环境温度',
                    '最大净功率',
                    '实际净功率',
                    '优化负荷',
                    '实际负荷',
                    '最优背压',
                    '实际背压',
                    '优化风机转速',
                    '实际风机转速',
                    '风机运行台数',
                    '优化风机功耗',
                    '实际风机功耗',
                    '调整后主汽流量',
                    '原始主汽流量',
                    '给煤量减少量',
                    '给煤节省成本 (元/分钟)',
                    '风机增加成本 (元/分钟)',
                    '经济效益 (元/分钟)',
                ])

                all_results_file = self.output_filepath.replace('.xlsx', '_所有结果.xlsx')
                try:
                    all_results_df.to_excel(all_results_file, index=False)
                    print(f"所有样本在所有风机转速下的完整结果已保存至：{all_results_file}")
                except Exception as e:
                    print(f"保存完整结果时出错: {e}")

                # 将最优结果保存到一个DataFrame中
                result_df1 = pd.DataFrame(self.result_list, columns=[
                    '样本编号',
                    '环境温度',
                    '最大净功率',
                    '实际净功率',
                    '优化负荷',
                    '实际负荷',
                    '最优背压',
                    '实际背压',
                    '优化风机转速',
                    '实际风机转速',
                    '风机运行台数',
                    '优化风机功耗',
                    '实际风机功耗',
                    '调整后主汽流量',
                    '原始主汽流量',
                    '给煤量减少量',
                    '给煤节省成本 (元/分钟)',
                    '风机增加成本 (元/分钟)',
                    '经济效益 (元/分钟)',
                ])

                print("优化结果维度为：", len(self.result_list))

                try:
                    result_df1.to_excel(self.output_filepath, index=False)
                    print(f"最优结果已保存至：{self.output_filepath}")
                except Exception as e:
                    raise IOError(f"保存 Excel 文件时出错: {e}")

    def plot_opt_results(self):
        df = pd.read_excel(self.output_filepath)

        # 获取样本编号作为 x 轴
        x = df['样本编号']

        # 图像输出路径
        output_dir = "./pics/优化图像"
        os.makedirs(output_dir, exist_ok=True)

        # 定义对比项：标题、优化列名、实际列名、保存文件名
        plot_items = [
            ("最大净功率对比", '最大净功率', '实际净功率', 'power_comparison'),
            ("负荷对比", '优化负荷', '实际负荷', 'load_comparison'),
            ("背压对比", '最优背压', '实际背压', 'backpressure_comparison'),
            ("风机转速对比", '优化风机转速', '实际风机转速', 'fan_speed_comparison'),
            ("风机功耗对比", '优化风机功耗', '实际风机功耗', 'fan_power_comparison'),
            ("主汽流量对比", '调整后主汽流量', '原始主汽流量', 'steam_flow_comparison'),
            ("成本对比", '给煤节省成本 (元/分钟)', '风机增加成本 (元/分钟)', 'cost_comparison'),
        ]

        # 绘图并保存
        for title, opt_col, real_col, filename in plot_items:
            plt.figure(figsize=(8, 5))
            plt.plot(x, df[opt_col], label='优化值', color='tab:blue', linestyle='-', marker='o')
            plt.plot(x, df[real_col], label='实际值', color='tab:red', linestyle='--', marker='x')
            plt.title(title)
            plt.xlabel("样本编号")
            plt.ylabel(opt_col.replace('(元/分钟)', '').strip())
            plt.legend()
            plt.grid(True)
            plt.tight_layout()
            plt.savefig(os.path.join(output_dir, f"{filename}_baseline.png"), dpi=300)
            plt.close()

    def Update_Turbine_Parameters(self, mech_in, p_real, Ne_k, Pc_k):
        """  使用梯度下降更新模型参数 ita_c 和 ita_t 以及流量系数列表
        :param mech_in: 当前输入
        :param p_real: 当前实际抽气压力(调节级+7个压力级)
        :param Ne_k: 当前时刻负荷
        :param Pc_k: 当前时刻背压
        :return:
        """
        # 使用当前参数运行模型
        self.Turbine.update_parameters(self.ita_c, self.ita_t, self.k_list)
        self.Condenser.update_parameters(self.kc, self.cf)

        states = self.Sim_Step(mech_in, Pc_k)
        Ne_sim = states[5]  # 仿真的机组负荷
        error_Ne = Ne_sim - Ne_k  # Ne误差

        # 数值梯度估计(有限差分法)
        delta_ita_c = self.ita_c / 1000
        delta_ita_t = self.ita_t / 1000

        # 计算 ita_c 的梯度
        ita_c_plus = self.ita_c + delta_ita_c
        self.Turbine.update_parameters(ita_c_plus, self.ita_t, self.k_list)
        states_plus_ita_c = self.Sim_Step(mech_in, Pc_k)
        Ne_plus = states_plus_ita_c[5]
        ita_c_minus = self.ita_c - delta_ita_c
        self.Turbine.update_parameters(ita_c_minus, self.ita_t, self.k_list)
        states_minus_ita_c = self.Sim_Step(mech_in, Pc_k)
        Ne_minus = states_minus_ita_c[5]
        gradient_ita_c = (Ne_plus - Ne_minus) / (2 * delta_ita_c)

        # 计算 ita_t 的梯度
        ita_t_plus = self.ita_t + delta_ita_t
        self.Turbine.update_parameters(self.ita_c, ita_t_plus, self.k_list)
        states_plus_ita_t = self.Sim_Step(mech_in, Pc_k)
        Ne_plus = states_plus_ita_t[5]
        ita_t_minus = self.ita_t - delta_ita_t
        self.Turbine.update_parameters(self.ita_c, ita_t_minus, self.k_list)
        states_minus_ita_t = self.Sim_Step(mech_in, Pc_k)
        Ne_minus = states_minus_ita_t[5]
        gradient_ita_t = (Ne_plus - Ne_minus) / (2 * delta_ita_t)

        # 更新参数
        self.ita_c -= self.learning_rate * error_Ne * gradient_ita_c
        self.ita_c = max(0.75, min(1.0, self.ita_c))  # 限制范围

        self.ita_t -= self.learning_rate * error_Ne * gradient_ita_t
        self.ita_t = max(0.75, min(1.0, self.ita_t))  # 限制范围

        Dt_sim, Pt_sim, Dc_sim, hc_sim, Tc_sim, Ne_sim, Pc_sim, Vf_sim, Dec_sim, Pc_next_sim = states
        Dj1_cur = Dt_sim
        # 计算压力级的实际抽汽流量
        D1_real, D2_real, D3_real, D4_real, D5_real, D6_real, D7_real = self.Di_calc(p_real[1:])
        # 实际的抽汽压力（调节级+7个压力级）
        P_real = p_real
        # 计算级间流量
        Dj2_real = Dj1_cur - D1_real
        Dj3_real = Dj2_real - D2_real
        Dj4_real = Dj3_real - D3_real
        Dj5_real = Dj4_real - D4_real
        Dj6_real = Dj5_real - D5_real
        Dj7_real = Dj6_real - D6_real
        D_real = np.asarray([[Dt_sim, Dj1_cur, Dj2_real, Dj3_real, Dj4_real, Dj5_real, Dj6_real, Dj7_real]])

        # 添加到缓冲区(滑动窗口设置为30min)
        if len(self.P_buf) < 30:
            self.P_buf.append(P_real)
            self.D_buf.append(D_real)
        else:
            # 移除最旧的数据
            self.P_buf.pop(0)
            self.D_buf.pop(0)
            self.P_buf.append(P_real)
            self.D_buf.append(D_real)
            # 将缓冲区转换为 NumPy 数组
            P_buf_array = np.vstack(self.P_buf)
            D_buf_array = np.vstack(self.D_buf)

            # 流量系数校正
            k0 = self.turbine_coff(D_buf_array[:, 0].reshape(1, -1), P_buf_array[:, 0].reshape(1, -1))
            k1 = self.turbine_coff(D_buf_array[:, 1].reshape(1, -1), P_buf_array[:, 1].reshape(1, -1))
            k2 = self.turbine_coff(D_buf_array[:, 2].reshape(1, -1), P_buf_array[:, 2].reshape(1, -1))
            k3 = self.turbine_coff(D_buf_array[:, 3].reshape(1, -1), P_buf_array[:, 3].reshape(1, -1))
            k4 = self.turbine_coff(D_buf_array[:, 4].reshape(1, -1), P_buf_array[:, 4].reshape(1, -1))
            k5 = self.turbine_coff(D_buf_array[:, 5].reshape(1, -1), P_buf_array[:, 5].reshape(1, -1))
            k6 = self.turbine_coff(D_buf_array[:, 6].reshape(1, -1), P_buf_array[:, 6].reshape(1, -1))
            k7 = self.turbine_coff(D_buf_array[:, 7].reshape(1, -1), P_buf_array[:, 7].reshape(1, -1))
            self.k_list = np.array([k0, k1, k2, k3, k4, k5, k6, k7]).reshape(-1)

    def Update_Condenser_Parameters(self, mech_pre, Pc_pre, Pc_cur):
        """
        :param mech_pre: 上一时刻的输入
        :param Pc_pre: 上一时刻的背压
        :param Pc_cur: 当前时刻的背压
        :return: 参数更新
        """
        # 使用当前参数运行模型
        self.Turbine.update_parameters(self.ita_c, self.ita_t, self.k_list)
        self.Condenser.update_parameters(self.kc, self.cf)

        mech_out = self.Sim_Step(mech_pre, Pc_pre)
        Pc_cur_pre = mech_out[-1] * 1000  # 背压(kPa)
        error_Pc = Pc_cur_pre - Pc_cur  # Pc误差

        # 计算 kc 的梯度
        delta_kc = self.kc / 1000
        kc_plus = self.kc + delta_kc
        self.Condenser.update_parameters(kc_plus, self.cf)
        states_plus_kc = self.Sim_Step(mech_pre, Pc_pre)
        Pc_plus = states_plus_kc[-1]
        kc_minus = self.kc - delta_kc
        self.Condenser.update_parameters(kc_minus, self.cf)
        states_minus_kc = self.Sim_Step(mech_pre, Pc_pre)
        Pc_minus = states_minus_kc[-1]
        gradient_kc = (Pc_plus - Pc_minus) / (2 * delta_kc)

        # 计算kec的梯度
        delta_cf = self.cf / 1000
        cf_plus = self.cf + delta_cf
        self.Condenser.update_parameters(self.kc, cf_plus)
        states_plus_cf = self.Sim_Step(mech_pre, Pc_pre)
        Pc_plus = states_plus_cf[-1]
        cf_minus = self.cf - delta_cf
        self.Condenser.update_parameters(self.kc, cf_minus)
        states_minus_cf = self.Sim_Step(mech_pre, Pc_pre)
        Pc_minus = states_minus_cf[-1]
        gradient_kec = (Pc_plus - Pc_minus) / (2 * delta_cf)

        # 更新参数
        self.kc -= self.learning_rate * error_Pc * gradient_kc
        self.kc = max(0.0, min(0.04, self.kc))  # 限制范围
        self.cf -= self.learning_rate * error_Pc * gradient_kec
        self.cf = max(0.7, min(1.0, self.cf))  # 限制范围


    @staticmethod
    # 根据抽汽压力计算实际抽汽流量(工程计算方法)
    def Di_calc(p_real):
        """ 计算每个压力级的抽汽流量 """
        Pa = 0.101  # 大气压
        EGJCQ = 3  # 抽汽损失百分比, 3 %
        pressures = p_real
        Di_real = []
        for p in pressures:
            ps_real = (p + Pa) * (1 - EGJCQ / 100) - Pa
            D_real = math.sqrt(p - ps_real) * (-0.4972 * p ** 3 + 7.227 * p ** 2 - 19.44 * p + 48.69)
            Di_real.append(D_real)
        return tuple(Di_real)

    @staticmethod
    # 汽轮机流量系数计算(滑动窗口最小二乘法)
    def turbine_coff(D, P):
        k = P @ np.linalg.pinv(D)
        return k

    def calculate(self):
        MAE_Ne = mean_absolute_error(self.Ne_data[:self.N], self.Ne_sim)
        MAPE_Ne = mean_absolute_percentage_error(self.Ne_data[:self.N], self.Ne_sim) * 100
        R2_Ne = r2_score(self.Ne_data[:self.N], self.Ne_sim)
        self.load_error = MAE_Ne
        print(f"机组负荷建模误差: MAE:{MAE_Ne:.2f}(MW) , MAPE:{MAPE_Ne:.2f}% ,R2:{R2_Ne:.2f}")
        MAE_Pc = mean_absolute_error(self.Pc_data[:self.N], self.Pc_sim)
        MAPE_Pc = mean_absolute_percentage_error(self.Pc_data[:self.N], self.Pc_sim) * 100
        R2_Pc = r2_score(self.Pc_data[:self.N], self.Pc_sim)
        self.pc_error = MAE_Pc
        print(f"背压建模误差: MAE:{MAE_Pc:.2f}(kPa) , MAPE:{MAPE_Pc:.2f} %, R2:{R2_Pc:.2f}")
        MAE_Pt = mean_absolute_error(self.Pt_data[:self.N], self.Pt_sim)
        MAPE_Pt = mean_absolute_percentage_error(self.Pt_data[:self.N], self.Pt_sim) * 100
        R2_Pt = r2_score(self.Pt_data[:self.N], self.Pt_sim)
        self.pt_error = MAE_Pt
        print(f"调节级压力建模误差: MAE:{MAE_Pt:.2f}(Mpa) , MAPE:{MAPE_Pt:.2f} %, R2:{R2_Pt:.2f}")
        MAE_Dec = mean_absolute_error(self.Dec_data[:self.N], self.Dec_sim)
        MAPE_Dec = mean_absolute_percentage_error(self.Dec_data[:self.N], self.Dec_sim) * 100
        R2_Dec = r2_score(self.Dec_data[:self.N], self.Dec_sim)
        self.dec_error = MAE_Dec
        print(f"凝结水建模误差: MAE:{MAE_Dec:.2f}(kg/s) , MAPE:{MAPE_Dec:.2f} %, R2:{R2_Dec:.2f}")

    def plot_sim_results(self):
        self.plt_scope(self.Ne_sim, '机组负荷(MW)', self.Ne_data[:self.N])
        self.plt_scope(self.Pc_sim, '背压(kPa)', self.Pc_data[:self.N])
        self.plt_scope(self.Pt_sim, '调节级压力(Mpa)', self.Pt_data[:self.N])
        self.plt_scope(self.Dec_sim, '凝结水流量(kg/s)', self.Dec_data[:self.N])

        self.plt_scope(self.Dc_sim, '末级排汽流量')
        self.plt_scope(self.hc_sim, '末级排汽焓值')
        self.plt_scope(self.kc_sim, '凝汽器还热系数')
        self.plt_scope(self.Rf_sim, '污垢热阻')
        self.plt_scope(self.ita_c_sim, '末级效率参数')
        self.plt_scope(self.ita_t_sim, '汽轮机效率参数')

    @staticmethod
    def plt_scope(y, title, y_true=None):
        # 图像输出路径
        output_dir = "./pics/仿真图像"
        os.makedirs(output_dir, exist_ok=True)
        size = 12
        plt.figure(figsize=(10.8, 3.2))
        plt.plot(y, color="r", linestyle="-", label="仿真值")
        if y_true is not None:
            plt.plot(y_true, color="b", linestyle="-.", label="真实值")
        plt.xlabel("t/(min)", fontsize=size)
        plt.ylabel(title, fontsize=size)
        plt.legend(loc="best", fontsize=size)
        plt.tick_params(axis='x', labelsize=size)
        plt.tick_params(axis='y', labelsize=size)
        plt.tight_layout()
        plt.grid()
        safe_title = title.replace("/", "_").replace("\\", "_").replace(":", "_")
        plt.savefig(os.path.join(output_dir, f"{safe_title}_baseline.png"), dpi=300)
        plt.close()
    def model_output(self):
        return {
            "计算机组负荷": self.Ne_sim,
            "计算背压": self.Pc_sim,
            "计算调节级压力": self.Pt_sim,
            "计算凝结水流量": self.Dec_sim,

            "末级排汽流量": self.Dc_sim,
            "末级排汽焓值": self.hc_sim,
            "凝汽器换热系数": self.kc_sim,
            "污垢热阻": self.Rf_sim,
            "末级效率参数": self.ita_t_sim,
            "汽轮机效率参数": self.ita_t_sim,

            "机组负荷建模误差": self.load_error,
            "背压建模误差": self.pc_error,
            "调节级压力建模误差": self.pt_error,
            "凝结水建模误差": self.dec_error,

            "最大净功率":self.opt_net_power,
            "优化负荷":self.opt_ne,
            "最优背压":self.opt_pc,
            "优化风机转速":self.opt_fan_speed,
            "优化风机功耗":self.opt_fan_power,
            "调整后主汽流量":self.opt_dst,
            "给煤量减少量":self.opt_coal,
            "给煤节省成本":self.opt_coal_cost,
            "风机增加成本":self.opt_fan_cost,
            "经济效益":self.opt_profit
        }


def get_data(ini_num, sim_num):
    path = r'.\data\8月数据-csv.csv'
    df = pd.read_csv(path, encoding='utf-8')

    data_NE = df['负荷'].values
    data_Pc = df['背压'].values
    data_Pst = df['主汽压力'].values
    data_Dst = df['主汽流量'].values
    data_Tst = df['主汽温度'].values
    data_Feed_water = df['主给水流量'].values
    data_Pre = df['再热蒸汽压力'].values
    data_Tre = df['再热蒸汽温度'].values
    data_T_Condwater = df['凝结水温度'].values
    data_D_Condwater = df['凝结水流量'].values
    data_Ta = df['环境温度'].values
    data_fan_speed = df['风机转速'].values
    data_fan_number = df['风机运行台数'].values
    data_Tc = df['排汽温度'].values
    data_Pt = df['调节级压力'].values
    data_Fan_P = df['风机功耗(转速)'].values
    data_Net_P = df['净功率'].values
    data_uB = df['给煤量'].values
    data_P_one = df['一段抽汽压力'].values
    data_P_two = df['二段抽汽压力'].values
    data_P_three = df['三段抽汽压力'].values
    data_P_four = df['四段抽汽压力'].values
    data_P_five = df['五段抽汽压力'].values
    data_P_six = df['六段抽汽压力'].values
    data_P_seven = df['七段抽汽压力'].values

    # 汽轮机模型输入
    Datas = np.column_stack((data_NE, data_Pc, data_Pst, data_Dst, data_Tst, data_Feed_water, data_Pre, data_Tre,
                            data_T_Condwater, data_D_Condwater, data_Ta, data_fan_speed, data_fan_number, data_Tc, data_Pt, data_Fan_P, data_Net_P, data_uB,
                            data_P_one, data_P_two, data_P_three, data_P_four, data_P_five, data_P_six, data_P_seven))
    Datas = Datas[ini_num:ini_num + sim_num, :]  # 9月份有效数据
    # Data = Data[:8600,:]              # 3月份有效数据
    # Tst, Pst, Dst, Pre, Tre, Ta, fan_speed, fan_number
    Uncertain_Datas = np.column_stack((data_NE, data_Pst, data_Dst, data_Tst, data_Pre, data_Tre,
                            data_Ta))
    # 虽然环境温度和AGC相关性不强（但用电和天气可能也存在关系），但是环境温度和AGC可能对主汽参数存在很大的影响，所以还是应该多输入多输出预测
    return Datas, Uncertain_Datas

if __name__ == '__main__':
    ini_num = 1440 * 7 + 12156
    sim_num = 2000
    Data, Uncertain_Datas = get_data(ini_num, sim_num)
    sim_model = Model_Simulation(Datas=Data, speed_change_limit=SPEED_CHANGE_LIMIT, speed_step=SPEED_STEP, output_filepath=OUTPUT_FILE_OPT, tolerance=TOLERANCE)
    sim_model.Run_sim_opt()
    sim_model.calculate()
    sim_model.plot_sim_results()
    sim_model.plot_opt_results()

    print('end')
